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Abstract 

The  goal  of  this  work  is  to  provide  a  method  for  choosing  joining  (e.g.,  bolt)  lo¬ 
cations  for  attaching  structural  reinforcements  onto  complex  structures.  The  join¬ 
ing  locations  affect  structural  performance  criteria  such  as  the  frequency  response 
and  the  static  compliance  of  the  modified  structure.  One  approach  to  finding  im¬ 
proved/optimal  joining  locations  is  to  place  the  joints  such  that  the  total  amount  of 
energy  input  into  the  structure  (from  external  forces)  is  lowered/minimized,  thus 
ensuring  that  the  performance  of  the  structure  is  least  affected  by  the  structural 
modifications.  However,  such  an  approach  does  not  account  for  the  stresses  in 
the  joints.  Therefore,  in  this  work,  the  amount  of  strain  energy  concentrated  in 
the  joints  is  also  considered.  The  cost  function  for  this  optimization  problem  is 
then  composed  of  two  energies.  These  energies  are  different  for  the  undamped 
and  damped  cases.  Herein,  the  focus  is  on  the  (more  realistic)  damped  case.  The 
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for  the  stresses  in  the  joints.  Therefore,  in  this  work,  the  amount  of  strain  energy  concentrated  in  the  joints 
is  also  considered.  The  cost  function  for  this  optimization  problem  is  then  composed  of  two  energies.  These 
energies  are  different  for  the  undamped  and  damped  cases.  Herein,  the  focus  is  on  the  (more  realistic) 
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cost  function  is  minimized  by  a  modified  optimality  criteria  method.  This  pro¬ 
cess  is  time  consuming  because  it  requires  the  calculation  of  sensitivities  of  the 
joint  strain  energy,  which  in  turn  requires  the  calculation  of  the  displacements  of 
all  candidate  joint  locations  by  using  the  system-level  mass  and  stiffness  matrices 
and  force  vector  (at  each  frequency  in  the  range  of  interest).  To  address  this  is¬ 
sue,  a  series  of  complex  algebraic  manipulations  and  approximations  are  used  to 
significantly  reduce  the  computational  cost.  In  addition,  for  the  case  where  struc¬ 
tural  and  geometrical  variations  are  necessary,  parametric  reduced-order  models 
are  used  to  compute  the  cost  function  with  further  significant  gains  in  computa¬ 
tional  speed.  Numerical  results  for  improved/optimal  joining  are  presented  for 
representative  complex  structures  with  structural  variabilities. 

Key  words:  joining  locations,  joint  strain  energy,  parametric  reduced-order 
models 

1.  INTRODUCTION 

Mechanical  structures  such  as  those  found  in  automobiles  and  airplanes  con¬ 
sist  of  multiple  components  which  are  assembled  using  joints  such  as  bolts,  welds, 
rivets,  etc.  The  locations  (assembly  points)  of  these  joints  affect  structural  perfor¬ 
mance  characteristics  such  as  the  static  compliance,  the  frequency  response,  and 
the  durability.  To  achieve  high  performance,  the  joining  locations  should  be  se¬ 
lected  by  a  systematic  approach  rather  than  an  experience-based  approach.  How¬ 
ever,  this  issue  can  be  quite  challenging  because  there  are  many  joints  and  even 
more  possible  joining  locations  for  large  scale  complex  structures.  The  number 
of  such  joints  can  be  as  many  as  several  thousand.  The  choice  for  joining  lo¬ 
cations  can  be  improved/optimized  by  topology  optimization  approaches  such  as 
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homogenization  techniques  [1]  and  density  methods  [2,  3,  4].  Homogenization 
techniques  compute  an  optimal  distribution  of  micro-structures  in  a  given  design 
domain.  Density  methods  compute  an  optimal  distribution  of  isotropic  materials, 
where  the  material  densities  are  design  variables.  Although  the  single-component 
topology  design  has  been  extensively  studied  during  the  past  two  decades  [5],  the 
amount  of  research  done  for  multiple-component  topology  optimization  is  rel¬ 
atively  small.  In  that  area  of  research,  Chirehdast  and  Jiang  [6]  extended  the 
concept  of  topology  optimization  to  the  design  of  spot-weld  and  adhesive  bond 
patterns.  A  year  later,  Jiang  and  Chirehdast  [7]  proposed  a  theoretical  framework 
to  determine  which  optimal  connection  points  minimize  the  static  compliance  of 
the  given  substructures.  To  solve  the  coupled  problem  of  component  topology  and 
joining  location  optimization,  Chickermane  and  Gea  [8]  considered  a  methodol¬ 
ogy  for  a  multiple-component  structure  as  a  whole,  in  which  the  optimal  topology 
and  the  joint  locations  were  computed  simultaneously.  More  recently,  Zhu  and 
Zhang  [9]  did  layout  optimization  of  structural  supports  using  a  topology  opti¬ 
mization  method  for  free  vibration  analyses.  All  these  previous  efforts  employed 
spring  elements  for  modeling  joints.  In  contrast,  Li  et  al.  [10]  proposed  a  fastener 
layout/topology  that  achieves  an  almost  uniform  stress  level  in  each  joint,  and 
adopted  evolutionary  structural  optimization  [11,  12,  13,  14,  15,  16]  to  provide 
an  alternative  optimization  strategy  to  traditional  gradient-based  topology  opti¬ 
mization  approaches.  In  the  context  of  these  past  efforts,  the  focus  of  this  work  is 
on  the  development  of  an  efficient  framework  for  determining  improved/optimal 
joining  locations  as  to  minimize  the  total  energy  input  into  the  structure  and  the 
strain  energy  in  the  joints  of  a  complex  structure  with  variability  using  a  density- 
based  method. 
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For  general  optimization  processes,  finite  element  models  (FEMs)  are  typ¬ 
ically  used  to  evaluate  the  cost  function.  However,  the  number  of  degrees  of 
freedom  (DOFs)  of  FEMs  of  complex  structures  are  prohibitively  large.  So,  con¬ 
ventional  FEMs  are  hard  to  employ  due  to  the  expensive  time  needed  for  each 
iteration.  To  reduce  the  computational  cost,  Craig-Bampton  component  mode 
synthesis  (CB-CMS)  was  employed  by  Ma  et  al.  [17]  in  a  multi-domain  topology 
optimization.  The  CB-CMS  method  is  one  of  the  most  well-established  methods 
for  constructing  reduced-order  models  (ROMs)  [18,  19,  20,  21,  22].  However,  if 
one  attempts  to  use  CB-CMS  techniques  when  parametric  changes  (such  as  thick¬ 
ness  and  geometrical  variations)  are  applied  during  the  design  or  exist  through 
damage,  the  ROMs  have  to  be  reconstructed.  This  reconstruction  requires  other 
analyses  in  addition  to  the  repetitive  calculation  of  the  cost  function.  This  is  com¬ 
putationally  expensive  and  requires  significant  effort  to  prepare  a  FEM  and  a  ROM 
for  each  reanalysis. 

These  challenges  are  addressed  in  this  work  as  follows.  First,  the  mean  com¬ 
pliance  for  the  dynamic  case  with  damping  is  derived,  and  the  strain  energy  in 
the  joints  is  added  to  the  cost  function.  Second,  a  novel  approach  to  calculate  the 
sensitivity  of  the  strain  energy  in  the  joints  efficiently  is  proposed.  Third,  the  cost 
function  and  its  sensitivity  are  computed  in  optimization  process  by  using  novel 
models  which  are  able  to  manage  structural  variabilities.  Recently,  design  oriented 
parametric  reduced-order  models  (PROMs)  have  been  developed  to  avoid  such 
prohibitively  expensive  reanalyses  of  complex  structures  [23,  24,  25,  26,  27,  28]. 
Here,  the  next-generation  PROMs  (NX-PROMs)  [28]  developed  by  the  authors 
are  employed  to  allow  complex  structures  to  be  divided  into  several  components 
when  determining  improved/optimal  joining  locations. 
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This  paper  is  organized  as  follows.  In  Section  2,  a  design  methodology  for  de¬ 
termining  improved/optimal  joining  locations  is  defined,  which  includes  models 
for  the  joints,  the  definition  of  the  associated  cost  function,  and  a  computationally 
efficient  method  to  determine  design  sensitivities  for  the  cost  function.  In  Sec¬ 
tion  3,  NX-PROMs  used  in  the  calculations  are  reviewed.  In  Section  4,  numerical 
simulations  are  used  to  demonstrate  the  proposed  approach  for  the  problem  of  at¬ 
taching  an  armor  plate  to  a  structure  with  a  V-shaped  bottom.  Finally,  conclusions 
are  summarized  in  Section  5 

2.  DESIGN  METHODOLOGY  FOR  OPTIMAL  JOINING  LOCATIONS 

In  single-component  topology  optimization,  the  primary  objective  is  to  obtain 
the  optimal  layout  of  the  structure.  When  multi-component  structures  are  con¬ 
sidered,  the  problem  is  extended  to  select  the  optimal  joining  locations  between 
components.  This  is  done  not  only  to  optimize  the  layout  of  each  of  the  subcom¬ 
ponents,  but  also  because  the  joining  locations  affect  the  structural  performance. 
Herein,  a  density-based  topology  optimization  technique  [2,  3,  4]  is  applied.  Roz- 
vany  et  al.  [29]  have  defined  this  method  as  a  modeling  technique  based  on  solid 
isotropic  material  with  penalization  (SIMP),  where  the  distribution  of  the  join¬ 
ing  stiffness  is  optimized  to  improve  the  static  or  dynamic  structural  performance 
of  the  entire  connected  structure.  The  SIMP  method  has  been  developed  to  re¬ 
place  the  size  and  orientation  variables  (of  the  holes  used  in  the  homogenization 
method  [1])  with  a  density  variable  (of  the  finite  elements)  in  the  design  domain. 
Herein,  the  idea  of  SIMP  is  employed  to  select  optimal  joining  locations  for  the 
entire  connected  structure.  To  improve/optimize  the  joining  locations,  the  stiff¬ 
nesses  of  the  joints  are  designed  using  density  functions.  Thus,  the  design  vari- 
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ables  are  the  densities  (or  stiffnesses)  of  the  joints.  These  densities  are  continuous 
variables  varying  between  0  and  1 .  A  location  where  the  joint  has  a  low  density 
(close  to  0)  is  not  effective/adequate  for  joining,  while  a  location  where  the  joint 
has  a  density  close  to  1  is  best  for  joining. 

2.1.  Design  Region  -  Models  for  Joints 

Being  one  of  the  controversial  tasks  in  structural  FE  analysis,  modeling  meth¬ 
ods  for  joints  have  been  extensively  studied.  Depending  on  the  required  accuracy 
and  complexity  of  the  problem  at  hand,  an  appropriate  modeling  strategy  can  be 
adopted  for  the  joints.  Several  techniques  for  modeling  joints  were  proposed  in 
the  literature  [30,  31,  32,  33,  34].  For  fatigue  analyses  based  on  local  stresses 
and  the  local  strength  of  the  material,  a  fine  mesh  of  the  structure  and  accurate 
joint  models  are  required.  For  noise  and  vibration  analysis  of  complex  structures, 
a  moderate  level  of  accuracy  and  complexity  is  required,  which  leads  to  simple 
models  for  the  joints.  For  a  design  optimization  problem,  and  especially  for  a  pre¬ 
liminary  design,  the  simple  flexible  bar  models  are  preferred  [6,  8,  32,  35]  because 
those  joint  models  can  be  easily  catered  toward  the  iterative  updating  employed  in 
design  optimization. 

In  this  work,  the  joints  are  modeled  as  three  rectilinear  springs.  Let  the  stiff¬ 
ness  associated  with  the  motion  of  one  of  the  two  ends  of  a  joint  (of  index  i)  be 
k,S:*  =  Diag  (  kXJ  kyA  kza  'j ,  where  kxa,  kyA.  and  kZJ  denote  the  stiffnesses 
of  the  spring  along  the  three  directions  of  a  local  Cartesian  reference  system  asso¬ 
ciated  with  joint  i.  Here,  Diag(v)  represents  a  diagonal  matrix  with  entries  given 
by  the  vector  v.  The  directional  stiffnesses  of  a  joint  are  often  related  to  each 
other.  In  this  work,  it  is  assumed  that  kX:i  =  kZA  =  otiki  and  kyA  =  ki.  This  is 
the  case  for  joints  such  as  bolts,  rivets,  spot  welds,  etc,  where  y  is  the  axis  of  the 
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joint  (e.g.,  the  axis  of  a  bolt).  Thus,  a  joint  is  modeled  as  having  6  DOFs  linked 
by  three  springs.  The  stiffness  matrix  for  the  ith  joint  can  thus  be  written  as 


Km  = 


where  ks  i  =  fc,Diag 


<%•  1  oij 


The  joints  (modeled  as  three  rectilinear  springs)  are  designed  using  density 
functions  in  the  SIMP  method.  According  to  the  SIMP  method,  the  design  ele¬ 
ments  are  written  using  the  densities  pi  as 


K 


b,i 


P 

Pi 


=  P?K 


(1) 


where  K/M  represents  the  joining  stiffness  matrix  for  the  ith  candidate  joining  lo¬ 
cation.  Thus,  Kfc  j  is  a  density-based  function.  Intermediate  values  of  pt  (0  <  pt  < 
1)  are  penalized  compared  to  values  of  0  or  1  by  the  use  of  the  penalty  exponent 
p.  This  exponent  is  typically  p  =  3  for  a  structural  optimization  problem  [36,  37]. 
Also,  for  simplicity,  we  assume  that  a \  has  the  same  value  for  all  joints.  Thus, 
only  one  density  variable  is  assigned  to  a  joint. 

The  system-level  joining  stiffness  matrix  K/,  is  given  by 


K(,  =  Bdiag 


KfM  K  bj2 


(2) 


where  Bdiag  denotes  a  block-diagonal  matrix,  and  g  is  the  number  of  candidate 
joining  locations.  Then,  the  system-level  governing  equation  for  the  structural 
dynamic  problem  with  (structural)  damping  7  is 


M0 


ub 

iir 


+  (1+77) 


Kb  0 
0  0 


ub 

ur 


(3) 
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where  j  =  \/—l,  M0  and  K0  are  the  system-level  mass  and  stiffness  matrices 
which  do  not  include  the  joining  stiffness  matrix  Kb.  Subscript  b  indicates  the 
candidate  joining  DOFs,  and  subscript  r  denotes  the  remaining  DOFs.  Note  that 
a  joining  mass  matrix  Mb  does  not  exist  because  massless  spring  elements  are 
used.  Based  on  Eq.  (3),  the  dynamic  response  of  all  DOFs,  ur  (remainder)  and  ub 
(joint),  are  obtained. 

2.2.  Formulation  of  the  Optimization  Problem 

The  approach  employed  here  is  based  on  an  energy  criterion  which  is  com¬ 
monly  used  in  structural  optimization  problems.  Two  energies  are  used.  The 
first  is  the  total  energy  input  into  the  structure  under  dynamic  loading.  This  total 
energy  input  is  equal  to  the  external  work  done  on  the  structure,  which  can  be  de¬ 
fined  in  function  of  the  mean  compliance  of  the  structure.  For  the  dynamic  case, 
the  total  work  done  on  the  structure  by  external  forces  is 

Re(^J  FTd\j^J  =  Re  ^  FT"^dt)  =  Re  Q  FTu d?j  ,  (4) 

where  F  =  f e]ujt  is  the  external  harmonic  forcing,  and  u  is  the  displacement  due 
to  the  harmonic  forcing.  The  phase  reference  for  the  calculation  is  chosen  such 
that  f  is  real.  For  structures  with  damping,  the  response  u  is  complex  and  can  be 
expressed  as 

u  =  (ur  +  juj)ejujt  so  that  ii  =  (jumR  —  uui)ejuit,  (5) 

where  subscripts  R  and  /  indicate  real  and  imaginary  parts.  From  Eq.  (5),  the  real 
valued  portion  of  the  velocity  is 

f?e(u)  =  —ujuIe:’ult  =  — tuu/costuf  —  uiuR  sin  ujt.  (6) 
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Substituting  Eq.  (6)  into  Eq.  (4)  and  using  the  fact  that  f  is  real,  one  obtains 


=  — 7rfTUj  =  — 7r(fTu)j. 


The  resulting  first  component  of  the  cost  function  c\  is  thus 

Ci  =  -(fTu)7,  (7) 

and  contains  the  strain  energy  in  the  entire  structure,  including  the  joints.  How¬ 
ever,  focusing  on  the  durability  of  the  joints,  the  strain  energy  in  the  joints  should 
be  taken  into  account.  Thus,  the  second  component  of  the  cost  function  is  based 
on  the  strain  energy  in  the  joints.  This  energy  can  be  expressed  as 

c2  =  ^ufKbufe,  (8) 

where  the  superscript  H  indicates  the  Hermitian  operator.  Then,  by  assembling 
the  two  components  Ci  and  c2  of  the  cost  function  from  Eq.  (7)  and  Eq.  (8),  the 
final  cost  function  for  this  optimization  problem  is 

\ 

C  =  W\C\  +  W2C2  =  ~W i  (f  1u)/  +  W2~ UfeK^Ufe, 

where  w\  and  w2  are  weighting  factors  to  control  the  relative  importance  of  overall 
structural  vibration  and  joint  durability. 

Naturally,  the  number  of  joints  to  be  distributed  in  the  design  domain  is  lim¬ 
ited.  Thus,  the  topology  optimization  problem  associated  with  the  joining  location 
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design  can  be  stated  as 


Minimize  :  c(p)  =  —w\  (fru)7  +  u^-UftK^Ufo, 

9  (9) 

Subject  to  :  g{p)  =  ^Pi  -  N  <0;  0<  pmin  <  pi  <  1, 

2=1 

where  N  denotes  the  total  number  of  joints  allowed  in  the  design,  g  is  the  total 
number  of  candidate  joint  locations,  and  pmin  is  a  sufficiently  small  lower  bound 
imposed  to  avoid  numerical  instabilities  (herein  pmin  =  0.001). 

To  solve  such  optimization  problems,  specific  methods  have  been  developed 
to  handle  a  large  number  of  design  variables  with  a  few  constraints.  Among  these 
techniques,  the  method  of  moving  asymptotes  (MMA)  [38,  39]  and  the  optimality 
criterion  (OC)  [1,  40]  methods  are  broadly  utilized  for  their  efficacy  and  gener¬ 
ality.  The  MMA  method  is  based  on  the  convex  approximation  method  with  the 
advanced  feature  of  setting  asymptotic  moving  limits  to  approximation  variables. 
The  OC  makes  use  of  the  well-known  Karush-Khun-Tucker  condition  to  satisfy  a 
set  of  criteria  related  to  the  behavior  of  the  structure.  Even  though  the  OC  method 
is  well-convergent  for  static  cases,  it  is  not  effective  for  the  dynamic  case.  Thus, 
herein  we  use  the  modified  optimality  criterion  (MOC)  method  [40],  which  is  also 
a  gradient-based  optimizer. 

2.3.  Sensitivities  of  the  Cost  Function 

The  design  domain  for  joining  is  modeled  with  the  density-based  three  rectilin¬ 
ear  springs,  each  having  a  design  variable  p^  (density),  as  in  Eq.  (1).  The  variable 
Pi  is  varied  between  0  and  1  using  the  MOC  method  to  select  improved/optimal 
joining  locations.  For  any  gradient-based  optimizer,  the  design  sensitivities  of  the 
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cost  function  and  of  the  constraints  with  respect  to  the  design  variables  are  re¬ 
quired.  For  an  efficient  calculation  of  the  design  sensitivities  for  the  dynamic  case 
discussed  here,  an  adjoint  variable  method  [41]  is  applied.  First,  we  consider  the 
design  sensitivities  of  ci  given  by  in  Eq.  (7).  The  derivative  of  c\  with  respect  to 
the  mth  design  variable  prn  is 
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The  direct  calculation  of 
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is  cumbersome,  so  an  adjoint  variable  A  is 


used.  To  obtain  A,  the  equilibrium  Eq.  (3)  is  differentiated  with  respect  to  the 
design  variable  prn  to  obtain 
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Multiplying  Eq.  (10)  by 
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Thus,  the  design  sensitivity  of  c\  to  prn  is  given  by 
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Second,  the  sensitivity  of  c 2  with  respect  to  prn  is  considered.  One  obtains 

dc2  (p) 
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(12) 


dpm.  dpm 

where  the  fact  that  Kb  is  a  real,  symmetric  matrix  was  used.  This  sensitivity 
requires  the  calculation  of  and  -f^.  The  first  term  can  be  easily  calculated 
because  it  has  a  simple  analytical  form.  Next,  from  Eq.  (10),  one  obtains 
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This  equation  could,  in  principle,  be  used  to  compute  J^-  once  is  calculated. 
However,  using  Eq.  (13)  requires  the  inverse  of  G  at  each  iteration.  This  matrix 
is  very  large  because  it  is  a  full-order,  system-level  matrix.  Also,  G  depends  on 
excitation  frequency  u,  so  this  inversion  has  to  be  done  at  each  frequency  in  the 
range  of  interest.  To  avoid  such  a  high  computational  effort,  we  propose  a  novel 
approach.  To  calculate  J^,  first,  Eq.  (3)  is  used  to  obtain 
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(14) 
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Then,  substituting  Eq.  (13)  into  Eq.  (14),  one  obtains 
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where  G0  =  (— tu2M0  +  (1  +  j'y)  K0).  The  quantity  G0G  1  in  Eq.  (15)  can  be 
written  as 


/ 


GnG  = 


I  +  (1  +  j'y) 


V 


Kb  0 
0  0 


-l 


Go1 


Then,  substituting  Eq.  (16)  into  Eq.  (15),  one  obtains 
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The  novel  approach  uses  the  assumption  that  the  values  of  the  (spring)  stiffnesses 
of  the  joints  are  much  smaller  than  the  values  of  the  stiffnesses  in  the  nominal 
structure.  Thus,  we  assume  that,  for  all  DOFs  of  indices  i\  and  i2. 
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Then,  the  inverse  term  in  Eq.  (17)  can  be  written  as 

K0  O' 

0  0 


1+  (1  +37) 


Go1 


I  -  (1+37) 


K0  0 
0  0 


Go1- 


Substituting  Eq.  (18)  into  Eq.  (17),  one  obtains 
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Of  course,  Eq.  (19)  can  be  used  to  obtain  J^-.  However,  Eq.  (19)  requires  the 
calculation  of  the  inverse  of  G0.  Inverting  G0  has  to  be  done  only  once  during 
the  iterations  because  G0  does  not  depend  on  Kfe.  However,  G0  does  depend  on 
the  excitation  frequency  u>,  which  means  that  Go  has  to  be  inverted  at  all  frequen¬ 
cies  in  the  range  of  interest.  Also,  Go  is  a  large  matrix  because  it  is  a  full-order, 
system-level  matrix.  Thus,  this  calculation  is  very  time  consuming.  A  new  ap¬ 
proach  to  address  this  issue  is  presented  next.  The  key  term  to  be  calculated  in 

i  r  §^ufe  i 

Eq.  (19)  is  G0  pr"  .  To  compute  this  term,  we  consider  first  that  a  unit 
0 

force  is  applied  at  the  mth  joining  location  and  to  the  ath  DOFs  (at  that  location). 
The  index  a  varies  from  1  to  the  number  L  of  DOFs  used  in  the  finite  elements 
which  contain  the  joining  node  m.  For  example,  L  —  6  for  shell-type  elements, 
while  L  —  3  for  brick-type  elements.  Then,  the  resulting  deformation  \kma  can 
be  expressed  as 

r  i T 

Vm,a  =  Gq  1  •  [  0  0  •  •  •  lm>0  0  0  J  , 

where  \&m)0  indicates  the  deformation  everywhere  in  the  system  due  to  the  unit 
force  applied  at  the  ath  DOFs  of  the  mth  joint  location. 
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Next,  we  express  the  difficult  term  as 
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Thus,  the  sensitivity  of  c2  in  Eq.  (12)  can  be  expressed  as 

dc2  (p)  1 
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2=1 


where 


Aj  =  pfu^K^  (1+  n) 
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a=  1 


Finally,  the  design  sensitivity  of  the  entire  cost  function  for  the  mth  design 
variable  is  obtained  using  Eqs.  (11)  and  (20)  as 

dc  (p) 


dpr 


=  -  Wl  (ppfn  1uJmKsufe>m)/ 

if  n 

+  W2-  -2  Yf,  Ai  +  P/ClufmKsut,, 


2=1 


Computing  the  sensitivities  of  the  entire  cost  function  based  on  this  formulation  is 
fast  especially  because  vectors  ^m,a  have  to  be  calculated  only  once.  They  remain 
unchanged  during  the  optimization  iterations. 
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3.  PARAMETRIC  REDUCED-ORDER  MODELS  (PROMs) 


The  cost  function  and  design  sensitivities  are  presented  in  Section  2  as  if  they 
are  calculated  based  on  full-order  finite  element  models.  However,  if  the  struc¬ 
ture  has  a  huge  number  of  DOFs,  the  turn-around  time  of  the  optimization  iter¬ 
ation  process  is  very  long.  This  issue  is  particularly  important  when  the  design 
involves  not  only  choosing  joining  locations,  but  also  modifications  of  various 
components  of  the  structure.  In  that  case,  the  full-order  model  of  the  modified 
components  changes,  which  requires  additional  computational  effort.  To  address 
this  issue,  a  new  modeling  approach  is  presented  next.  This  approach  is  based 
on  next-generation  parametric  reduced-order  models  (NX-PROMs)  [28]  used  to¬ 
gether  with  the  well-known  fixed-interface  Craig-Bampton  component  mode  syn¬ 
thesis  (CB-CMS)  [19].  A  review  of  this  approach  is  provided  next. 


3.1.  Modeling  Approach 

CB-CMS  [19]  is  used  to  model  only  the  substructures  which  do  not  have  any 
structural  variability.  This  modeling  approach  is  used  because  it  is  very  simple  and 
computationally  stable.  To  apply  CB-CMS,  the  complex  structure  of  interest  is 
divided  into  several  substructures,  and  their  DOFs  are  partitioned  into  internal  and 
interface  DOFs.  The  interface  DOFs  for  a  substructure  (of  index  q)  are  projected 
onto  the  generalized  coordinates  by  using  static  constraint  modes  .  The  internal 
DOFs  are  projected  onto  fixed-interface  normal  modes  .  Then,  the  size  of  the 
mass  and  stiffness  matrices  and  the  force  vector  for  substructure  q  is  significantly 
reduced  as  follows 
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where  the  superscript  C  indicates  generalized  interface  DOFs  (i.e.,  constraint  par¬ 
titions).  These  DOFs  are  used  to  assemble  substructural  matrices  and  obtain 
system-level  reduced  matrices.  The  superscript  N  indicates  generalized  internal 
DOFs.  These  DOFs  are  used  to  reduce  the  number  of  internal  DOFs. 

The  substructures  which  can  have  variability  are  modeled  using  NX-PROMs  [28]. 
One  important  advantage  of  NX-PROMs  is  that  the  finite  element  mesh  of  the 
nominal  structure  does  not  need  to  be  modified  although  several  substructures 
may  have  variability.  That  is  because  the  mass  and  stiffness  matrices  of  these  sub¬ 
structures  are  parameterized.  The  NX-PROM  approach  resembles  the  CB-CMS 
approach.  However,  the  transformation  matrix  for  NX-PROMs  is  constructed  for 
all  values  of  the  variable  parameters  in  the  parameter  space  of  each  component 
with  variability.  In  contrast,  components  with  no  parameter  variability  do  not 
need  a  parameterization,  so  they  are  modeled  by  CB-CMS.  By  applying  the  NX- 
PROM  approach  to  the  Ith  substructure  with  variation  A p  in  one  of  its  parameters, 
the  mass  and  stiffness  matrices  and  force  vector  are  obtained  as 
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3.2.  Geometric  Compatibility  Conditions 

The  complete,  reduced-order  component-level  equations  of  motion  for  each 
component  /  of  the  entire  set  of  n  components  can  be  expressed  as 

MfOMqfOM  +  KfOMqfOM  =  FfOM,  (21) 

where  the  superscript  ROM  indicates  that  either  CB-CMS  or  NX-PROM  was 
used,  with  qi  being  the  generalized  coordinates  (/  =  !,•••,  n). 
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The  constraint  partitions  (indicated  by  superscript  C )  of  component-level  ma¬ 
trices  retain  the  physical  meaning  of  the  interface  DOFs.  This  means  that  the 
geometric  compatibility  conditions  at  the  interfaces  with  no  joints  can  be  applied 
directly  to  construct  the  system-level  matrices.  Consider,  for  example  that  an  in¬ 
terface  with  no  joints  exists  between  components  l  and  d  (d  —  1,  •  •  •  ,  n,  d  ^  l). 
Then, 


qp  =  qp> 


(22) 


where  qp  and  qp  are  the  generalized  coordinates  for  the  constraint  partitions  that 
correspond  to  the  interface  between  substructures  /  and  d.  Of  course,  there  is  no 
compatibility  condition  to  be  enforced  for  two  components  which  do  not  have  a 
common  interface.  Equation  (22)  is  used  to  transform  the  matrices  in  Eq.  (21)  in 
a  manner  similar  to  the  assembly  process  in  all  finite  element  modeling  methods. 
Then,  the  system-level  equation  of  motion  which  does  not  include  the  joints  is 
given  by 


1 \t ROM  "ROM  I  tsROM  ROM  p ROM 

•LVJ- sys  '  *^sys  ^i.sys  r  sys 


(23) 


Equation  (23)  is  obtained  after  all  geometric  compatibility  condition  have  been 
enforced,  except  for  the  conditions  present  at  the  interfaces  with  joints.  To  tackle 
the  joints,  the  (remaining)  constraint  partitions  corresponding  to  the  joints  are 
repartitioned  in  two  pieces.  These  pieces  are  indicated  by  superscript  G)  and  C2. 
The  C\  portion  corresponds  the  DOFs  of  one  end  of  all  joints  and  the  C2  portion 
corresponds  to  the  DOFs  of  the  other  end  of  all  joints.  Thus,  the  matrices  in 
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Eq.  (23)  can  be  expressed  as 
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where  superscript  C  represents  the  constraint  partition  (for  all  components)  that 
corresponds  to  the  joints. 

Next,  the  joints  (three  rectilinear  springs)  are  applied  to  connect  the  DOFs  of 
Ci  to  those  of  C2.  First,  the  joining  stiffness  matrix  in  Eq.  (2)  and  the  physical 
coordinates  of  all  DOFs  of  all  joints  are  partitioned  similar  to  Cj  and  C2  to  obtain 
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The  C\  and  C2  partitions  are  the  same  in  Eqs.  (24)  and  (25).  Thus, 


(25) 


qCl  =  qf1  and  q°2  =  qf2. 


(26) 


Ultimately,  Eq.  (26)  is  used  to  obtain  the  final  system-level  equation  of  motion 

with  joints  expressed  as 
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where 


Msys 


F 


sys 


KCl  +  K^1  K^lC2  KCliV 

kc2C!  kc2  +  kc2  kc2n 

J^NCi  J^jVC'2  J£iV 


and 


F 


sys- 


Note  that  all  the  design  parameters  are  contained  in  the  joints.  Thus,  the  join¬ 
ing  optimization  has  to  reevaluate  only  the  joints,  and  does  not  require  a  reevalu¬ 
ation  of  all  components.  Moreover,  all  components  (except  for  the  joints)  are  re¬ 
duced  only  once,  at  the  initial  construction  of  NX-PROMs,  before  iteration.  Thus, 
the  joining  design  can  be  very  efficient  by  using  NX-PROMs  with  the  proper  ma¬ 
trix  partitioning.  By  using  the  system  matrices  in  Eq.  (27)  (based  on  NX-PROMs), 
the  tum-around  time  of  the  iteration  process  is  much  shorter  than  by  using  FEMs. 
Additionally,  variations  in  any  substructure  (where  NX-PROM  is  used)  can  be 
handled  efficiently  in  the  new  optimization  processes. 


4.  NUMERICAL  EXAMPLE:  V-SHAPED  BOX  STRUCTURE  WITH  THICK¬ 
NESS  VARIATIONS 

To  demonstrate  the  improved/optimal  joining,  a  structure  with  a  V-shaped  bot¬ 
tom  is  considered,  as  shown  in  Figure  1.  The  focus  of  this  example  is  to  find 
joining  locations  where  to  attach  an  armor  plate  to  the  structure.  Figure  1  shows 
all  substructures  and  their  number.  The  marked  regions  are  candidate  joining  lo¬ 
cations.  Harmonic  loads  are  assumed  to  act  on  substructure  5  shown  in  Figure  3. 
Substructure  6  and  the  armor  plate  have  thickness  variations.  Table  1  shows  two 
cases  of  thickness  variation  of  each  substructure.  NX-PROMs  are  constructed  by 
parameterizing  the  thickness  of  substructure  6  and  of  the  armor  plate.  CB-CMS 


20 


is  applied  to  all  remaining  substructures  because  they  do  not  have  structural  vari¬ 
ations. 

As  an  initial  guess,  all  the  candidate  joining  nodes  on  substructure  4  and  on  the 
armor  plate  are  connected  by  three  rectilinear  springs  as  shown  in  Figure  2.  For  all 
joints,  the  maximum  allowable  stiffness  of  the  spring  in  the  (main)  y  direction  is 
ky  =  k0  =  500  kN/m,  and  the  other  directional  stiffnesses  are  kx  =  kz  =  0.5  k0. 
The  total  number  of  candidate  joints  is  54,  and  the  final  desired  number  of  joints 
g  is  10  or  11.  Note  that  the  four  edge  nodes  highlighted  in  Figure  2  are  not 
considered  candidate  joining  locations.  The  fact  that  joints  are  present  at  those 
four  locations  is  considered  to  be  known. 

The  optimization  starts  with  an  initial  guess,  i.e.  a  given  set  of  feasible  design 
values  {pi  =  0.185  when  N  =  10,  and  pt  =  0.204  when  N  =  11).  The  structural 
damping  is  7  =  0.03,  and  the  weighting  factors  in  Eq.  (9)  are  w\  =  0.5  and 
w2  =  0.5.  With  the  given  initial  guess,  the  excitation  frequency  was  fixed  at 
30  Hz  for  the  nominal  structure.  Figure  4  shows  the  results  of  the  optimization 
for  the  30  Hz  excitation  and  N  =  10.  Figures  5,  6  and  7  show  the  results  of  the 
optimization  for  a  100  Hz  excitation  for  the  nominal  structure  and  for  cases  1  and 
2  of  thickness  variations.  For  the  nominal  structure  under  the  100  Hz  excitation, 
N  —  11.  For  cases  1  and  2  of  thickness  variation,  N  —  10. 

The  optimal  joining  locations  described  in  Figure  5  are  selected  differently  for 
each  case,  even  though  the  thickness  variations  are  not  very  large.  However,  the 
cost  function  is  minimized  for  all  3  cases,  as  shown  in  Figure  6.  Figure  7  shows 
the  changes  in  natural  frequencies  at  each  iteration.  Note  that  the  various  choices 
made  during  the  iterations  affect  significantly  the  dynamic  response  because  some 
of  the  natural  frequencies  of  the  overall  structure  change  significantly. 


21 


5.  CONCLUSIONS 


Several  challenges  of  current  methods  for  determining  improved/optimal  join¬ 
ing  locations  have  been  addressed.  First,  the  mean  compliance  for  the  dynamic 
case  with  damping  was  derived,  and  the  strain  energy  in  the  joints  was  added  to 
the  cost  function.  Second,  a  novel  approach  to  calculate  efficiently  the  sensitiv¬ 
ity  of  the  strain  energy  in  the  joints  was  proposed.  Third,  the  cost  function  and 
its  sensitivities  were  computed  in  the  optimization  process  by  using  novel  next- 
generation  parametric  reduced-order  models  to  improve  computational  efficiency 
and  to  manage  structural  variabilities  in  several  substructures. 

The  approach  to  select  improved/optimal  joining  locations  uses  a  density- 
based  topology  optimization  method  which  employs  solid  isotropic  material  with 
penalization  (SIMP)  modeling.  Based  on  SIMP  modeling,  a  three  rectilinear 
springs  (with  density)  is  used  to  model  each  joint.  Also,  a  reliable  cost  func¬ 
tion  has  been  developed.  It  includes  the  energy  input  into  the  structure  and  the 
strain  energy  in  all  joints.  By  penalizing  the  density  of  the  springs  between  0 
and  1,  the  cost  function  is  minimized  while  satisfying  a  constraint  which  enforces 
an  upper  limit  for  the  number  of  joints  in  the  design  domain.  To  solve  this  opti¬ 
mization  problem,  the  modified  optimality  criterion  method  has  been  applied.  To 
demonstrate  the  methodology,  the  problem  of  attaching  armor  to  a  structure  with 
a  V-shaped  bottom  has  been  considered.  By  applying  the  proposed  methodology, 
improved/optimal  joining  locations  have  been  selected. 
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Table  1 :  Thickness  variations  for  substructure  5  and  for  the  armor  plate 


Substructure 

Nominal 

Case  1 

Case  2 

6 

6  mm 

7.5  mm 

8.5  mm 

Armor  plate 

10  mm 

10.5  mm 

11.1  mm 
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Candidate  joining  locations 


Figure  1 :  Structure  with  a  V-shaped  bottom;  indices  of  the  substructures  are  shown 
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/ 

Three  directional  springs 


Figure  2:  Springs  modeling  joints  between  substructure  4  and  the  armor  plate 
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excite  symmetric  modes 


excite  asymmetric  modes 


Figure  3:  Dynamic  loads  applied  to  substructure  5  to  excite  symmetric  and  asymmetric  modes  of 
the  entire  structure 
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(a) 


(b) 


(c) 

Figure  4:  (a)  10  optimal  joining  locations,  (b)  convergence  history,  and  (c)  natural  frequency 
variations  for  a  30  Hz  excitation  of  the  nominal  structure 
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(C) 

Figure  5:  Optimal  joining  locations  for  a  100  Hz  excitation  for  (a)  nominal  structure,  (b)  case  1  (c) 
case  2  of  thickness  variation 
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(a) 


(b) 


(c) 

Figure  6:  Convergence  history  for  a  100  Hz  excitation  for  (a)  nominal  structure,  (b)  case  1  (c)  case 
2  of  thickness  variation 
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(C) 

Figure  7:  Natural  frequency  variation  for  a  100  Hz  excitation  for  (a)  nominal  structure,  (b)  case  1 
(c)  case  2  of  thickness  variation 
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